#ifndef __QGCM_PROFILE_READERS_SRC__
#define __QGCM_PROFILE_READERS_SRC__

#include <limits>
#include "readers.h"


REAL fourier_transform_coefficient = numeric_limits<REAL>::quiet_NaN();
REAL fourier_transform_coefficient_squared = numeric_limits<REAL>::quiet_NaN();

void set_fourier_transform_coefficient(const REAL & val){
	fourier_transform_coefficient = val;
	fourier_transform_coefficient_squared = val * val;
};

void computecrossprodmagnitude(REAL& rval, const REAL & k, const REAL & p, const REAL & q){
	REAL k2 = k*k; REAL p2 = p*p; REAL q2 = q*q;
	REAL f1 = (k2+p2+q2);
	REAL f2 = 2.0*(k2*k2+p2*p2+q2*q2);
	REAL ifactor = f1*f1-f2;
	rval = 0.5*sqrt(ifactor);
}

void computelawsinescoeff(REAL & rval, const REAL & k,const REAL & p, const REAL & q){
	computecrossprodmagnitude(rval, k, p, q);
	rval = rval/(k*p*q);
};


qgcomplex compute_asymKPQintegrand(const REAL & k, const REAL & p, const REAL & q, const INDEX & np, const INDEX & mk){
	// due to problems with the asin being ambiguous when computing angles that are obuse in the triangle (i.e. theta > iI/2)
	// we have switched this function to use law of cosines to compute the 2 interior angles needed
	REAL k2 = k*k;
	REAL p2 = p*p;
	REAL q2 = q*q;
	REAL cosbeta = (k2+q2-p2)/(2.0*k*q);
	REAL cosgamma= (k2+p2-q2)/(2.0*k*p);
	REAL gamma = acos(cosgamma);
	REAL beta  = acos(cosbeta);
	const REAL npf = static_cast<REAL>(np);
	const REAL mkf = static_cast<REAL>(mk);
	REAL expval = ( npf*gamma - (mkf - npf)*beta );
	//REAL expval2= ( npf*gamma + (mkf - npf)*beta );
	// return p*q*qgcomplex(cos(expval), sin(expval));
	//return p*q*qgcomplex(0.5*(cos(expval)+cos(expval2)), 0.0); // ignore the sin part as a test for now
	return qgcomplex(cos(expval), 0.0)*p*q; // ignore the sin part as a test for now
};

bool computeP_def(const rank3 & findex){
	const model_class & mo = getModel();
	static INDEX atmos_layer_1_index = mo.get_atmosphere_layer_1_field_number();
	static INDEX atmos_layer_2_index = mo.get_atmosphere_layer_2_field_number();
	static INDEX ocean_layer_1_index = mo.get_ocean_layer_1_field_number();
	static INDEX ocean_layer_2_index = mo.get_ocean_layer_2_field_number();
	static INDEX mixed_layer_index   = mo.get_mixed_layer_field_number();
	static INDEX topo_layer_index    = mo.get_topo_layer_field_number();
	if(findex==rank3(atmos_layer_1_index,atmos_layer_1_index,atmos_layer_1_index)) return true;
	if(findex==rank3(atmos_layer_2_index,atmos_layer_2_index,atmos_layer_2_index)) return true;
	if(findex==rank3(ocean_layer_1_index,ocean_layer_1_index,ocean_layer_1_index)) return true;
	if(findex==rank3(ocean_layer_2_index,ocean_layer_2_index,ocean_layer_2_index)) return true;
	if(findex==rank3(ocean_layer_2_index,topo_layer_index,ocean_layer_2_index)) return true;
	if(findex==rank3(mixed_layer_index,mixed_layer_index,ocean_layer_1_index)) return true;
	return false; 
};

//-------------------------------------------------------------------------------------------------
// REAL computeQ
//-------------------------------------------------------------------------------------------------
REAL computeQ_def(const model_class & mo, const REAL & k, const rank2 & qindex){
	// called the first tiem this function is called
	static vector<REAL> fcoeffs = mo.get_global_prof_by_name(gdb_prof_str_FCOEFFS).getData();
	static INDEX atmos_layer_1_index = mo.get_atmosphere_layer_1_field_number();
	static INDEX atmos_layer_2_index = mo.get_atmosphere_layer_2_field_number();
	static INDEX ocean_layer_1_index = mo.get_ocean_layer_1_field_number();
	static INDEX ocean_layer_2_index = mo.get_ocean_layer_2_field_number();
	static INDEX mixed_layer_index   = mo.get_mixed_layer_field_number();
	static INDEX topo_layer_index    = mo.get_topo_layer_field_number();

	const REAL k2 = k*k;
	const REAL nb = - fourier_transform_coefficient_squared*k2;
	// #define nb (-39.4784176*k2) // n stands for negative b
	if( qindex == rank2(atmos_layer_1_index, atmos_layer_1_index) ) return nb - fcoeffs[0];
	if( qindex == rank2(atmos_layer_1_index, atmos_layer_2_index) ) return fcoeffs[0];
	if( qindex == rank2(atmos_layer_2_index, atmos_layer_1_index) ) return fcoeffs[1];
	if( qindex == rank2(atmos_layer_2_index, atmos_layer_2_index) ) return nb - fcoeffs[1];
	if( qindex == rank2(ocean_layer_1_index, ocean_layer_1_index) ) return nb + fcoeffs[2];
	if( qindex == rank2(ocean_layer_1_index, ocean_layer_2_index) ) return -fcoeffs[2];
	if( qindex == rank2(ocean_layer_2_index, ocean_layer_1_index) ) return -fcoeffs[3];
	if( qindex == rank2(ocean_layer_2_index, ocean_layer_2_index) ) return nb + fcoeffs[3];
	if( qindex == rank2(mixed_layer_index  , mixed_layer_index) ) return 1;
	if( qindex == rank2(topo_layer_index   , topo_layer_index) ) return 1;
	return 0;   
#undef nb
}

REAL computeA_wmap (const rank6 & pos, const int nfs, const map_container<REAL, 3> * Q,
		const map_container<REAL, 3> * QI, 
		bool (*computeP) (const rank3 & ) ){
	double rval =0.0;
	const INDEX & i = pos[3];
	const INDEX & j = pos[4];
	const INDEX & k = pos[5];
	const INDEX & ki = pos[0];
	const INDEX & pi = pos[1];
	const INDEX & qi = pos[2];
	// Q, QI maps are have all indices starting at 0
	// A_{ijk} = sum_{m} Q^{-1}_{im} \sum_{n} ( P_{mnk} Q_{nj} + P_{mnj} Q_{nk} )
	for(INDEX m = 1; m <= nfs; m++){ // sum over m !!
		const double & qival = QI->getvalue( rank3( ki, i-1, m-1 ) );
		double isum = 0.0;
		for(INDEX n = 1; n <= nfs; n++){
			if( computeP( rank3( m, n, k ) ) ){
				const REAL & rqval = Q->getvalue( rank3(pi, n-1, j-1) );
				const double qval = static_cast<double>(rqval);
				isum += qval ;
			} // symmetrization is done automatically since all elements are sampled by the static computation 
			if( computeP( rank3( m, n, j ) ) ){
				const REAL & rqval = Q->getvalue( rank3(qi, n-1, k-1) );
				const double qval = static_cast<double>(rqval);
				isum -= qval;
			}
		}
		rval += qival*isum;
	}
	rval *= -(fourier_transform_coefficient_squared)*0.5;
	//rval *= (-39.4784176)*0.5; // 0.5*(2*PI*i)^2 // comes from the fourier transform of the Jacobian 
	return static_cast<REAL>(rval);
};
//-------------------------------------------------------------------------------------------------
#endif
